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The reliability and precision of numerically solving stochastic non-Markovian equations by stan- 
dard numerical codes, more specifically, with the fourth-order Runge-Kutta routine for solving 
differential equations, is gauged by comparing the results obtained from analytical solutions for the 
equations. The results for different prescriptions for transforming the non-Markovian equations in 
a system of Markovian ones are compared so to check the reliability of the numerical method. 
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I. INTRODUCTION 

Stochastic equations of motion and theirs generahzations are extensively used in different contexts, e.g. in classical 
statistical mechanics to study systems with dissipation and noise, to determine how order parameters equilibrate, in 
critical phenomena dynamics, among many other problems [ij. 

The simplest case of a stochastic equation of motion is the phenomenological Langevin equation describing e.g. 
the Brownian motion. It consists of local (Markovian) dissipation and noise terms, which are related to each other 
through the classical Fluctuation-Dissipation theorem. Memory effects are then neglected. However, in real systems 
scattering events responsible to dissipation and fluctuations proceed through finite time intervals, which, consequently, 
result in finite memory effects for these quantities. The typical equations of motion describing real physical systems 
are then expected to be nonlocal (i.e., non-Markovian) ones with memory effects. An equation of this type is given 
by a generalized Langevin equation (GLE) of motion of the form (for a general reviews, see e.g. Ref. Q), 

4>{t) + I dt'K{t - t')^{t') + v'icb) = m , (1.1) 



where is a variable of the system (for example the coordinate of a particle) in interaction with a thermal bath (the 
dot means derivative with respect to time), V{(l)) is a potential term that is a function of </> (with prime denoting 
derivative with respect to (p), K{t — t') is the dissipation kernel and the noise term ^{t) is a Gaussian fluctuation 
with zero mean but colored, i.e., with two-point correlation satisfying the generalized classical Fluctuation-Dissipation 
relation at temperature T (throughout this work we consider the Boltzmann constant equal to one), 

mm)^TK{t-t'). (1.2) 

Models with equations of the form of Eq. (|l.ip are also known as Caldeira-Leggett type of models Q . Applications 
making use of equations of the form of Eq. (|l.ip must be able to properly deal with the memory kernel. In some restrict 
cases, like in the most common forms of non-Markovian kernels used in the literature, e.g. when the dissipation kernel 
K{t — t'), describes an Ornstein-Uhlenbcck (OU) process [1] or in the case of the Exponential Damped Harmonic 
(EDH) kernel , it is possible to reduce the non-Markovian equation into a set of Markovian ones with white noise 
properties. In other generic cases, however, this may not always be possible (see Ref. [6| for a recent review on the 
different colored noises and associated equations used in the literature). In these and any other cases dealing for 
example, with nonlinear equations, we must resort to numerical methods. Though there are some specific numerical 
methods that may be applicable for general cases 0, we still would like to be able to solve equations like Eq. p.ip 
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through standard methods, which are less numerically expensive than other alternatives. Analytically, if the equation 
is linear, then it is, in principle, possible to solve equations like Eq. (jl.ip through a Laplace transform, since the 
dissipation integral term in Eq. (|l.ip is just in the form of a convolution. Even so, in these cases we can only look at 
averaged (over the noise) quantities. In addition, if it is nonlinear, i.e., when the potential term V{(j)) is a polynomial 
form of order larger than two in then we must resort to numerical methods in order to solve the differential 
stochastic equation of motion. In the cases it can be reduced to a set of local differential equations with white noise, 
like in the OU and EDH cases mentioned above, the most natural way for solving the system of differential equations 
would be e.g., through a standard Runge-Kutta method, which is both easy to implement and usually produces results 
with good accuracy. However, since the standard Runge-Kutta method is basically a deterministic algorithm, when 
applied to a stochastic differential equation with white noise, it becomes not well defined, since the white noise term 
has infinite variance and, therefore, cannot be generated. In order to deal with this problem, stochastic Runge-Kutta 
routines have been proposed @ . An issue related to applying such techniques to solve systems of differential equations 
is that not all equations may be stochastic, but only a sub-set of them. In cases like these, those algorithms may not 
be suitable or appropriate. 

In this work our objective is to investigate both analytically and numerically the solutions of non-Markovian linear 
equations of the form of Eq. (jl.ip and then compare their results. In doing so, we are able to gauge the reliability and 
precision of numerically solving the stochastic non-Markovian equations by standard numerical codes, e.g., Runge- 
Kutta codes for solving differential equations. For this, we use the most common forms for the dissipation kernel 
K{t-t') given by the OU and EDH forms. 

The paper is organized as follows. In Section 2, we briefly describe the solution of the linear GLE through Laplace 
transform. In Section 3, we show the transformation of the GLEs of motion with OU and EDH kernels in systems of 
Markovian time differential equations with white noise. In Section 4, we show the comparison of the results obtained 
for the analytical solutions for the linear equations with those obtained numerically within our prescription to make 
then local, from a standard fourth-order Runge-Kutta code. Finally, in Section 5 we present our conclusions and final 
comments about the precision of numerically solving the equations with standard numerical methods. 

II. THE GENERALIZED LANGEVIN EQUATION: LINEAR REGIME 

Since non-Markovian equations like Eq. have nonlocal kernel terms in the form of a convolution, they become 
suitable to be solved by Laplace transform. If we write the potential V{(j)) in the form 



m = — + (2.1) 

where is a parameter of the potential and we have separated the interaction term (non-quadratic) V/(0) from the 
quadratic one. By neglecting interaction terms in Eq. (|2.ip . the GLE Eq. p.ip takes the linear form, 

'(i>{t) + m^(t>{t) / dt'K{t - t')<j){t') = ^{t) . (2.2) 
Jo 

By making use of the Laplace transform for (j){t), 



dt exp(— st)0(t) 



(2.3) 



and from the convolution theorem applied to the non-Markovian dissipation term in Eq. 
that the solution for the linear GLE can be written in the Laplace transform form as 



we can easily obtain 



m 



K{s) m 



sK{s) 



sK{s) 



(2.4) 



where K{s) and ^(s) are the Laplace transforms of the dissipation kernel K{t — t') and the noise £,{t), respectively. 
The solution for (/)(t) is obtained from the inverse transform of Eq. 
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Jo 



(2.5) 



where 



0(0)+ s + i^(s) 0(0) 



sK{s) 



(2.6) 



and 



g(t - t') = 



sX(s) 



(2.7) 



The exphcit solution for (p{t) is difRcuh to give analyticaUy because of the noise term on the right hand side of 
Eq. (j2.5p . but since the noise is Gaussian, (^) =0, we obtain that its average is simply given by 



(2.8) 



and once the kernel K{t — t') is given, it is easily computed through Eq. (|2.6p . either numerically or algebraically. 
We here have used the MAPLE software to numerically evaluate for (p{t). It should be noted that in the OU and 
EDH cases the explicit forms for the solutions can be obtained by MAPLE, but they are too complicated and long 
solutions, so we refrain ourselves here to write them down explicitly. 

It is also convenient to calculate {(fP{t)). Remembering that the non-Markovian noise ^{t) satisfies Eq. (II. 2p . 
{^{t)^{t')} = TK{t - t'), we then also obtain that 



'-{t)) =ip\t)+T f dt"g{t-t") [ dt'g{t-t')K{t' ~t") 
Jo Jo 



(2.9) 



III. THE OU AND EDH KERNEL CASES 

Let us now describe the two cases of dissipation/noise kernels we are interested in studying here, as mentioned in 
the introduction, the OU and EDH cases, which are also the most common forms of non-Markovian kernels used in 
the literature. We will describe the equations with these two types of kernels separately. 

A. The GLE with OU kernel 

Many studies considering the influence of Gaussian colored noise on nonlinear physical systems are usually made 
considering the OU noise, with two-point correlation satisfying 



i^oumouit')) = TKou{t-t') 



(3.1) 



with 



-7(t-*') 



(3.2) 



where 7 gives the inverse of the time scale for the kernel memory and Q is the overall magnitude of the dissipation, 
which also gives the magnitude of the dissipation in the local (Markovian) limit ^ , 



POO 

Q= dt'K{t - t') 
Jo 



(3.3) 
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It can be easily shown that the OU noise can be generated by the stationary part of the solution of the following 
differential equation: 



iou{t)^~^[iou{t)-V^C\ , (3.4) 
where C in Eq. p.4p is a white Gaussian noise satisfying 

(CW) = , 

(am')) = 6{t-t'), (3.5) 

If we now define the new variable Wqu (t) j given by 

t 



Wou {t) = - dt'Kou {t - t') 4>{t') , (3.6) 
Jo 

it can be shown that Wqu satisfies the equation of motion 

Wou (t) = -jWou (t) ~ Kou (0) {t') . (3.7) 

Using Eqs. p.4p and p.7p . we can transform the integro-diffcrcntial Eq. p.ip with OU kernel into a fourth-order 
dynamical system of local first-order differential equations given by 

= ?/ , 

y - -V'{4>)+iou + Wou , 

Wou = ~-fWou~Kou{Q)4> , 

iou = -7 [Coc/ - y2Tg C] • (3.8) 

B. The GLE with EDH kernel 

Another case of non-Markovian equation of interest in the literature is the one with a EDH kernel, whose noise 
term, ^^(i), satisfies 

{iH{t)) = , 

{(.Hmnit')) = TKnit-t') , (3.9) 

with kernel Kh {t — t') given by 

ifij(t-0 = e"^'*"*'^^{cos[r!i(t-i')]+j^sin[^^i(t -<')]} , (3.10) 

where Q and 7 have the same meaning as in the OU case and fif = fig — 7^ > 0. It can be easily shown that the 
noise can be generated by the following differential equation 's]: 

'init) + 2-fiH{t) + nl^Hit) = nly/2TQC{t) , (3.11) 



where (^(t) is a white Gaussian noise with the same properties as given in the OU noise case, expressed by Eq. p.Sp . 

Similarly as in the OU case, by defining a new variable Wh analogous to Eq. (|3.6[) and after analogous algebra 
leading to the system of equations (j3.8|) . we obtain that the GLE with EDH kernel can be written in terms of a 
sixth-order dynamical system of local first-order differential equations given by 
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= y , 

y = -V'{cp) + Wh+£.h , 

Wh = u-2-fWH - KHiO)y , 

in = z , 



where we have also defined a new function u{t) as 



(3.12) 



dt' 



dK„(t-t') 
dp 



-2-fKH{t-t') 



d(j){t') 
dt' 



(3.13) 



C. Alternative Prescription 



A different prescription than the one leading to the system of differential equations (j3.8|) and (|3.12|) is, instead of 
defining the function W{t) like in Eq. p.6p (and similarly with the EDH kernel case), would be to define it as 



W{t) 



dt'K {t -t')(P{t') +^{t) 



(3.14) 



whose only difference with the previous prescription is the addition of the noise term to the equation. This prescription 
is used with some frequency in the literature, e.g., like in 

In terms of Eq. (|3.14p . the system of differential equations, Eq. (|3.8p . for the OU case then changes to^ 



y 

Wou 

while Eq. (|3.12p for the EDH case changes to 



-r (0) + Wou , 

7\/2rQ C - iWou 



Kou (0) y 



(3.15) 



y 

Wh 

u 

z 



y , 

-V'{(b) + Wh , 

u-2j{WH-CH)-KH{0)y + z , 
-nliWn - Ch) + KH{0)y - 2-fKHiO)y , 
z , 

-2-fz-nl^H + nl^/2fQC ■ 



(3.16) 



The two systems of differential equations, Eqs. (|3.15p and p.l6p . are equivalent to the previous two, Eqs. (|3.8p and 
p.l6p . and they produce identical results. The main difference between the two prescriptions being the fact that 
the latter requires some extra care in its numerical implementation. Since the function W, as defined by Eq. (I3.14p . 
involves the noise function £,{t), its initial condition must be carefully set in terms of stationary solution of the noise 
differential equations, Eqs. (|3.4p and p. lip , for the OU and EDH cases, respectively, otherwise the resulting dynamics 
between the two different prescriptions for deriving the two sets of differential equations will lead to different results. 



^ Note that in the OU case, with the choice Eq. 113.141 1. the noise ^ou decouples from the system and only the white noise contribution 
remains. 
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IV. COMPARING ANALYTICAL AND NUMERICAL SOLUTIONS 

Now we show our analytical results for {(j){t)) and {(/^{t)) obtained using Laplace transformation and compare with 
our numerical results obtained by solving the systems of equations derived before from the two prescriptions used in 
the previous section, i.e., Eqs. p.Sp and (|3.12p and Eqs. (|3.15p and (|3.16p . As commented at the end of the last 
section, the two prescriptions lead to identical results. For the results shown below we used the first prescription, 
leading to the systems of differential equations, Eqs. p.Sp and (|3.12p . 



(a) (b) 




121 ^ 1 ^ 1 ^ 1 . 1 ^ 1 

10 20 30 40 50 
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Figure 1: Time evolution for ip(t) in the OU case: (a) for 7 — 0.5, (b) for 7 — 1.0 and (c) for 7 = 5.0. The other parameters are 
taken as m = 1.0, Q = 1.0 and T = 1.0. 

In Fig. [1] we plot side by side our results for (p{t) obtained from the analytical expression Eq. (|2.6[) and those 
obtained numerically by solving the system of first-order differential equations, Eq. (13. 8p . for the OU case. In Fig. [2] 
the same is done for the system Eq. p.l2p for the EDH case, including the analytical solution for this same case. The 
system of differential equations (|3.8p and (|3.12p are solved by a standard fourth-order Runge-Kutta algorithm with 
time stepsize of At = 0.01. The number of realizations over the noise used in both OU and EDH cases was 300, 000. 
In all our simulations we have used the initial conditions (t>(0) = 1 and 4'{0) = 0. 

From both Figs. [1] and [5] we see an excellent agreement between the analytical and numerical results obtained from 
the standard Runge-Kutta code. For completeness it is also useful to compare the results for {(/P), defined analytically 
by Eq. (|2.9p . with the numerical results for this same quantity. This is done in Fig. [3] for the OU (left panel) and 
EDH (right panel) cases, respectively. 

From Fig. [3] we again see an excellent agreement between the results obtained for ((/)^) analytically and numerically. 
How good is this agreement between analytical and numerical results in both cases can also be better assessed by 
defining the difference between them, i.e., 

A(f> 

— V-'analytic V^numeric 
A(j)^ = ((^^)analytic - (0^)numcric ■ (4.1) 
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(a) (b) 




I • \ • \ • \ • \ • I _n 4 1 • I • I • I I L 

10 20 30 40 50 10 20 30 40 



t t 



(b) 




Figure 2: Time evolution for ifl{t) in the EDH case: (a) for 7 — 0.1, (b) for 7 = 0.3 and (c) for 7 = 0.5. The other parameters are 
taken as flo — 




Figure 3: The time evolution for {(j>^{t)) in the OU case (left panel) and EDH case (right panel). The parameters used are: 7 = 0.5, 
Qo = 1.0, m = 1.0, Q = 1.0 and T = 1.0. 



The results for the differences and are shown in Figs. |3]and[5l respectively, for the OU and EDH cases. 

We note from the results shown in Figs. S] and [H] that the differences are allways smaller than about 10~^ and 
oscillates around zero in a noisy way. In fact we have checked that most of this difference is purely due to noise and 
can be decreased by increasing the number of realizations over the noise. This certifies that the solution from the 
standard fourth-order Runge-Kutta algorithm for the system of differential equations (|3.8p and (|3.12p is reproducing 
quite well the analytical results, despite the initial non-deterministic character of the GLE. The agreement is seen 
both at short times, where the memory effects dominate, but also at long times, where it becomes sub-dominant and 
where the local approximation with dissipation Eq. (|3.3p can better represent the dynamics [l^, [ll| . 
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Figure 4: The difference A(/) in the OU case (left panel) and EDH case (right panel). The parameters used are: 7 = 0.5, Qo = 1-0, 

m=1.0,Q = 
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Figure 5: The difference Acfy^ in the OU case (left panel) and EDH case (right panel). The parameters used are: 7 = 0.5, Qo = 1-0, 
m = 1.0, Q = 1.0 and T = 1.0. 

We think that the overall error observed between the analytical and numerical results can probably be made even 
smaller with an improved code, like for example by using a stochastic Runge-Kutta one [8]. However, we were not 
able to fully verify it for the particular cases of stochastic differential equations studied here, since its use would mean 
solving all equations in (19) and (23), or (26) and (27), the same form, which seems not appropriate, since they are not 
all stochastic. We hope to better discriminate this problem in a future work, where a variation of the Runge-Kutta 
code is in test to be used in situations like these. 



CONCLUSIONS 



In this work we have studied the reliability of using standard numerical codes to solve generalized Langevin equa- 
tions. In this study we have used a standard fourth-order Runge-Kutta routine to solve for the generated system of 
local first-order differential equations. We have shown that the solution for the linear equation of motion obtained 
from the use of a Laplace transform, when contrasted with the numerical solution arc in very good agreement and the 
use of these standard numerical methods can lead to a reliable description of the stochastic dynamics, independent 
of the form of the prescription used to transform the original generalized Langevin equation in a system of local 
differential equations. 

We have studied the two most used cases of dissipation/noise kernels, the OU and EDH cases. We have observed 
that the results between the analytical and numerical ones agree with each other with very good numerical precision. 
We expect that the numerical precision can be made even better by using variations of a stochastic Rung-Kutta code, 
appropriately tailored to deal with systems of local differential equations like the ones we here have studied. Work in 
this direction is in progress and we hope to report on the results in a future publication. 

Though we have studied in this work only the OU and EDH non-Markovian cases, our results are expected to be 
of importance for the practical study of other different generalized Langevin equations through numerical methods. 
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which is the case in most situations, like when including nonlinear effects in the equations. In those cases, once a 
proper prescription is used to transform these equations, our results show that standard numerical methods for solving 
differential equations can be applied to obtain reliable results for the dynamics at both short and long times. 
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